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Algebraic-based nonstandard time-stepping schemes 


Loic MICHEL 


Abstract 

In this preliminary worlQ we present nonstandard time-stepping strate¬ 
gies to solve differential equations based on the algebraic estimation method 
applied to the estimation of time-derivative, which provides interesting 
properties of ’’internal” filtering. We consider firstly a classical finite dif¬ 
ference method, like the explicit Euler method for which we study the 
possibility of using the algebraic estimation of derivatives instead of the 
usual finite difference to compute the numerical derivation. Then, we in¬ 
vestigate how to use the algebraic estimation of derivatives in order to 
improve the slope predictions in RK-based schemes. 
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1 Introduction 


Ordinary differential equations and stiff differential equations DEI have been studied 
extensively and successful methods have been proposed (e.g. lailllSlEllTlIHlElEQ]), 
including the solver routines ’ode45’ and ’ode23s ’ dn na [i3]. Recent alternative 
methods propose for instance to use specific series, like power series [H] or Haar 
wavelets da [IS] to describe the solution of ODEs; the purpose is to substitute the 
standard finite difference by the numerical properties of the series. In this paper, we 
propose to extend the hnite difference technique by an estimation of the derivatives 
using the algebraic estimation approach. Introduced in ra. the algebraic estimation 
method [IHl US] has been widely applied to many different problems of estimation that 
occur in dynamical systems. Some of these applications aim e.g. to reconstruct the 
states of dynamical systems [201 HH], and a computational toolbox has been released 
to help processing efficiently the algebraic estimation of derivatives for particular 
problems [2T] . The algebraic estimation of derivatives could be seen as similar to the 
differentiation by integration technique, investigated especially by |22| and [23], and 
probabilistic Kriging-based method [2T| . 

Our proposed method belongs a priori to the class of nonstandard finite-difference 
(NSFD) methods [25l[26],|27l|28l[29] , described as ’’powerful numerical methods that 
preserve signihcant properties of exact solutions of the corresponding differential equa¬ 
tions” (an interesting survey can be found in [30]); explicit rules to ’’design” NSFD 
schemes have been proposed in [25] [26]. Among the derivations that have been pro¬ 
posed (e.g. [3I11321I331E1]), one emphasizes the contribution of [3l] that describes 
a nonstandard hnite-difference scheme for fractional systems, that uses a discrete 
version of the Caputo fractional derivation. 

In this paper, we attempt to design simple NSFD schemes, that are derived from 
the algebraic estimation technique in order to evaluate the discrete derivatives of first 
order ODEs considering an ’’Euler framework” and a ”RK framework”. In an ’’Euler- 
like scheme”, the algebraic estimation of the derivatives is mainly used to build a 
multi-step scheme, where the evaluation of the derivative depends also on the past 
values of the solution. In a RK-like scheme, the estimation of the local slope, that 
is usually ’’measured” as an average of several local slopes, is performed using the 
properties of hltering of the algebraic technique. 

The paper is structured as follow. Section 2 described the proposed methods, 
including a brief review of the algebraic estimation method. Some concluding remarks 
can be found in Section 3. 
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2 Outline of the method 

Consider an ordinary differential equation (ODE) such as: 


= f{y{t),u{t)), t e [0, tf], y{0) = Vo. (1) 

The quantities u and y represent respectively the input and the solution of ([^. The 
corresponding usual discrete explicit Euler scheme reads: 


dyjt) 

dt 


k-\-l 


Vk+I - Vk 

h 


~ f{yk,Uk), 


keN, y{0) = yo 


( 2 ) 


where k is the sampled time, h is the time-step, and yk, yk+i are respectively the 
solution of (|^ at the discrete instants tj. and The sampled are supposed equally 
distributed i.e. tk — tk-i is constant for all k. 


We identihed two major strategies to build time-stepping schemes, based on the 
rules that help designing NSED schemes and the algebraic estimation technique. The 
”Euler-like”(or ’’multi-step”) strategy follows the principle of the Euler method, and 
the ”RK-like” strategy is based on the general Runge-Kutta method [12]. 


Algebraic estimator of derivatives 

Using the previous notations, consider a function g{x) and dehne the algebraic esti¬ 
mator of derivative of order n such as: 


dg{x) 


dx 


T^{9)n = 


k-\-l 


®o5'fc+i + o^igk -l- ■ • • -l- angk-r 

(p{h) 


(3) 


0(h) = h + as h —)■ 0 where h = Xk — Xk-i is constant 

The coefficients aQ, Oi, • • • , are real coefficients. These coefficients and the function 
0(h) are dehned by ([W) in the proposition (2.1) as calculating steps of the proposed 
time-stepping schem^ We have: 


ao = T 

Oj = 2(T — 2hi), i E {I ■ ■ - 7] — 1} 

an = {T - 2hr]) (4) 

1 3Kh 3Kh 


^The algebraic estimator of derivative is defined with the highest degree of g that is equal to 
fc -I- 1. Depending on the application, this degree can be decreased to k. 
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2.1 Euler-like Algebraic-NSFD scheme 

2.1.1 Definition 

The proposed ”Euler-like” Algebraic-NSFD scheme aims to extend the finite difference 
in Q by the algebraic estimator 'D(/)„. First, one proposes a ’’symbolic” nonstandard 
scheme, that is of the form: 


i.e. 


^^(/)n ~ f{yk,Uk) 
(^OVk+l + OilVk + ■ ■ ■ + anUk-n+l 


(5) 


m 


e I/O = 2/(0). 


A forward scheme can be easily deduced from the proposed general scheme (|^: 


O^lVk + ■ • • + OtnUk—n+l 0(^) rf \ ; _ 1 ^ 1 *+ 

yk+i = -^- f{yk,Uk), fc e N+, I/O = 1/(0). (6) 

ao Oo 

As presented in j34] regarding the Caputo fractional derivative, due to the nonlocal 
nature of the algebraic-based derivative operator, the discrete representation of the 
derivative must take into account a part of the past history of the solution. The 
number n of the involved sampled solutions yk, - ■ ■ , yk-n defines a window that char¬ 
acterizes the ’’precision” of the derivative estimation|^ 


Then, in the following proposition, we formalize the proposed Fuler-like nonstan¬ 
dard time-stepping scheme based on the algebraic estimation framework. 

Proposition 2.1. Consider the following nonstandard numerical scheme associated 
to the ODE that verifies the general scheme 0; 

3h f 1 

K—<^ao + J22yk-r,+j+i{T-2jh)j^f{yk,Uk), keN*^, yo = y{0) 
with T = rjh > 0] rj ^ and ag = yk-ri+iT + yk+i{T — 2r/h) 


where k is the sampled time, h is the time-step, K is a real constant, and T is a 
multiple of h that characterizes the ’’low filtering” property of the algebraic derivative 
This scheme is called Euler-Algebraic-NSFD scheme, or simply 
E-A-NSFD scheme with a window T. 

•^This estimation window implies that the n initial sampled solutions are obviously not known at 
the beginning of the algorithm. To initialize the estimation window, one may consider e.g. using 
the classical finite difference scheme. 


(see §i in 2.1.2 
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Proof. Hypothesis We consider solving the ODE Q, for which the solution y{t) is 
assumed, in the time domain, to be locally represented by a linear function of the 
time i.e.: 


y{t) = ao + ait (8) 

where, in particular, the coefficient ai is calculated from the algebraic estimation 
technique. The lowest degree of ”time-£ltering” is considered. 


Step 1 - Algebraic Derivation The general technique to perform the estimation of 
derivatives using the algebraic estimation strategy is mainly described in mm- 
It allows estimating the coefficient ai in ([^. We consider only the hrst algebraic 
derivative of / according to the dehnition of the initial-value problem Q. Transform 
hrst ([^ in the Laplace domain, then take the derivative d/ds'. 

y{s) = — + ^ sy{s) =ao + — (9) 

s s 


y{s) + s 


d^ 
d s 



( 10 ) 


Step 2 - Back to the time domain Multiply hrst (10) by s 

s~^y{s) + 

ds 


( 11 ) 


then using the Cauchy formula]^ applied to each term of ( |TT| ) : 

r Ip - rr-'i-rVMr)dr ( 12 ) 

we deduc^the expression of ai that corresponds to the derivative of y through 


ai = 


(p- 1)! Jt- 


\LP ~ l p - ry-^Tx{T)dT 




(13) 


(p -h 1)! Jt-i 

Since we consider the lowest degree of ”time-£ltering”, we take p = 2. 


"^Some elements of proof can be found in [3S] (pp.22-23). 

®The expression of ai corresponds to the direct application of (12), whose integral is evaluated 
over a moving windows T. 
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Therefore, (13) reads: 


6 /•* , 


2r)|/(r)(i r 


(14) 


Step 3 - Integral expression of the numerical scheme Replacing in (j^ the time- 
derivative operator by (14) gives an integral expression of the proposed E-A-NSFD 
scheme: 


- -2T)y(T)<iT K f{y{t),u(t)) 


(15) 


Step 4 - Finally, the discretized version... The usual Trapezoidal scheme allows in¬ 
tegrating numerically and thus providing a discretized expression of (15). We have: 


Oi = —K 


3h 


ri-l 1 

Vk-n+iT yk+i{T - 2r]h) + ^ 2yk-n+j+i{T - 2jh) \ 


(16) 


which, as a result, gives the proposed NSFD scheme ([^ when substituted in ([^ ac¬ 
cording to the rule ^2 of the design methodology of NSFD scheme^ [30]. Depending 
on the value of T, the total gain of |ai| could be different, that is why we multiply oi 
by a positive constant K that should be estimated at the initial instant^ Note that 
if = 1, then T = h which cancel the summation term 2yk- “ 2jh). 


Step 5 - In addition, identihcation with the general scheme The coefficients at, i & 
{1 ■ ■ ■ n} can be identihed to the general form ([^. We have: 


ao = T 

OLi = 2(T — 2hi), i G {1 ■ ■ - 7] — 1} 

an = {T- 2hy) (17) 

1 3Kh 3Kh 


As an example of the Prop. 2.1, the simplest version of the proposed E-A-NSFD 
is a particular case where T = h that verihes ([^. 


® Since the complete finite difference scheme in ([^ is a priori replaced by a weighted sum of the 
discrete values yk, we deduce that the rule #2, regarding the design of NSFD scheme [30], is not 
completely verified (since only the ”h” term has to be changed). A future work should investigate 
this ambiguity. 

^In simulation, it is observed that iF —> 1 when T increases. 
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To illustrate how the E-A-NSFD scheme is written, consider the simple ODE: 


dyjt) 

dt 


byit) +u{t), 


2 /( 0 ) = 1/0 


( 18 ) 


Considering rj = 3] h 
(18) reads: 


10 ® (T 3 is a real constant), the E-A-NSFD scheme of 


T 3 ( 3 Mfc -2 + 2 |/fc_i - 2yk - 3yk+i) ^ 5y{t) -F 50u{t) 


(19) 


Considering rj = 5] h = 10 ^ (Ts is a real constant), the E-A-NSFD scheme of 


(18) reads: 

d' 5 ( 52 /fc -4 + Qyk -3 + 22 /fc -2 - 2i/fc_i - 6yk - 5yk+i) ^ 5y{t) + 50u{t) (20) 


As stated in Prop. | 2 . 2 | , the distribution of the signs in these two cases is equally 
distributed. 


2.1.2 Some properties 


This strategy could been seen as a receding horizon approach because the interval of 
integration T corresponds to the window from which the estimation of the derivative 
is performed. Therefore, the constant T should be chosen small enough in order to 
evaluate the estimate within an acceptable short delay, but large enough, in order to 
ensure a good low pass hltering [ 21 ]. 

Since the whole set of coefficients Oj does not depend of the discrete solution y^, 
they need only to be evaluated at the initialization step of the scheme. 

The following simple proposition gives an interesting property of the general 
scheme (|^ concerning the distribution of the signs of the set of a*. 


Proposition 2.2. The distribution of the signs (over the whole set of the coefficients 
ai, i = 0---r]) verifies: sign(Q;o, • ■ ■ , aLW 2 -ij) = -sign(Q;L7?/2j, • ■ ■ , a??)- 


Proof. From Prop, 
number of a 

*< LIJ. 


2.1 


we have a* = 7 (T — 2hi), i = 0 ■ ■ ■ 7 , 7 G N*’ 
that are e.g. positive, check the inequality (T — 2hi) > 


To get the 
This gives 

□ 


The E-A-NSFD scheme may become unstable for high 7 ]^ 

®The window T of the algebraic estimation may be seen as an ’’internal” time-delay inside the ODE 
(§ that may ’’deform” the value of the time-derivative estimation. Further works may investigate 
the conditions of stability according to the length of the window T and the coefficient K. 
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2.1.3 Toward possible generalizations 

One may consider a possible generalization of the scheme Q using a more general 
expression 0(h) instead of the constant time-step h. The scheme can be rewritten: 


with 0(h) 


met) 

J'2, 


77-1 f 

2yk-^+j+i{T - 2j0) I 




h + 0{h‘^), as h —)■ 0. 


( 21 ) 


2.2 RK-like Algebraic-NSFD scheme 

The family of RK methods is described generally as an Euler scheme for which mul¬ 
tiple estimations of the local slopes around f{yk) are computed in such manner that 
the solution increment yk+i is based upon a weighted average of these multiple esti¬ 
mations. The general RK scheme reads: 


Vk+I = yk + hY^ biQi (22) 

i=l 

where gi, q 2 -, - ■ ■ ,qs are the estimations of the different slopes through /. The number 
s characterizes the choice of the method in the RK family. 


2.2.1 Definition 


The proposed ”RK-like” Algebraic-NSFD scheme aims to substitute the weighted sum 
X)i=i biQi in (22) by the algebraic estimator V{f)n. 

First, one proposes a ’’symbolic” nonstandard scheme, that is of the form: 


yk+i = yk + (j){h)f{yk + 0(h)[/(i/fc) -F 0(h)T>(/)„]) i.e. 

0 V) ^ 4>{h)[f{yk) + 4>{h){aof{yk + ho) H- K ^nfiyk + hn))]) (23) 

k e K*+, yo = y{0) 

where the coefficients ho) hi, • • • , h„ are real coefficients such as for all i G [0, n], h* G 
[hmm, hmax]; «o ''' Q^n are real coefficients and 0(h) = h -|- Oih?), as h —)■ 0. 

Then, in the following proposition, we formalize the proposed RK-like nonstandard 
time-stepping scheme based on the algebraic estimation framework. 
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Proposition 2.3. Consider the following nonstandard numerical scheme associated 
to the ODE that verifies the general scheme 


Vk+i - Vk 

fiQi) 


= / I 2/fc + ^ 


fiVk) + 


mh^ 

T3 


ri-l 

<^o + 2/(2/fc + dj){T - 2jh) 

j=i 


with T = r]h>0-,r] e and ao = f{yk + do)T + f{yk + h^)(T - 2r]h), 

k e N*+, yo = y{0) 


(24) 


where k is the sampled time, h is the time-step, K is a real constant, and T is a 
multiple of h that characterizes the ’’low filtering” property of the algebraic derivative 
(see §i in 2.1.2). This scheme is called RK-Algebraic-NSFD scheme, or simply RK- 
A-NSFD scheme with a window T. 


Proof. Hypothesis We consider solving the ODE (Q, for which the fnnction / i 
assnmed, to be locally represented by a linear fnnction i.e.: 


IS 


f{x) = ao + aix (25) 

where, in particnlar, the coefficient ai is calcnlated from the algebraic estimation 
techniqne. The lowest degree of ’’time-filtering” is considered. The pnrpose is to 
evalnate the local valne of Oi for x G [xk-i, Xk], for all k > 1. 


Step 1 - From the discretized version of the derivative estimation From Prop, 
the expression of Oi in the discrete domain is given by: 


2.1 


ai = -K 


3h 


rj-l f 

f{xk-n)T + fixk)iT - 2r]h) -h 51 ‘^f{xk-ri+j)iT - 2jh) | 


(26) 


The set of valnes Xk-rj, Xk-rj+i, ■ ■ ■ ,Xk is, according to ([^, theoretically a regular grid, 
but since the interval [xk-i, Xk], A: > 1 is very small, we can consider evaluating p + l 
values of f{yk + Sj) where for all i G [0, -|- 1], hj G [Smin, dmax]- In particular, the 6i 

coefficients can be random numbers. We have: 


ai = —K 


3h 


T7-1 

f{xk + 6o)T + f{xk + 5,)(T - 2r]h) + V{xk + - 2jh) 

i=i 
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Step 2 - Definition of an RK2 scheme with algebraic estimation of / To build the 
proposed NSFD scheme, we start from the standard RK2 scheme. 

Considering the evolution of the solution y between the steps yk and yk+i, for any 
k > 1. The hrst step of the RK2 scheme, called the prediction, is to evaluate the 
solution yk+i /2 in the middle of the step i.e.: 


yk+l/2 — Uk + 

that gives, via (|^, the derivative of y at the middle of the time-step: 


dt 


k+l /2 


— fiyk+1/2) — f (yk + • 


(28) 


(29) 


Then, the second step, called the correction, is to evaluate the solution at the end of 
the step yk taking into account the derivative ^ 


fc+1/2 


i.e.: 


yu+i = yk + hf{yk+i/2)- (30) 

Instead of considering a single evaluation of the derivative at the middle of the time- 
step (or multiple evaluations of the derivatives like e.g. in the RK4 case), the proposed 
modihcation of the RK2 strategy is to estimate multiple local values of the slopes via 
/ between the steps yk and yk+i using the algebraic estimator V[f)n- First, one 
generalizes the ’’prediction” step (28) that is rewritten: 


yk+1/2 = l/fc + 2 < f^Vk) >h (31) 

where we denote < f{yk) >h the averaged value of f{y) between the steps k and k + 1 
(in other words, over the time-step h, regarding the notations)|^ 

Before rewriting the ’’correction” step, one needs the formal dehnition of the deriva¬ 
tive of a function g{x) taken at the point x = a: 


(h 

dx 


= lini 

h^O 


g{a + h)- g{a) _ g{a + h) - g{a) 
h ^ h 


for h very small (32) 


that allows dehning the quantity < f{yk) >h, by substituting (27) in (32): 


< f{yk) f{yk) + 


df{y) 


dy 


h = f{yk) + hai 


(33) 


Finally, from (31) and (33), the correction step is thus rewrittenj^ 

^Regarding the RK4 algorithm (with no time dependence), we have: < f{yk) >h= \{f{yk) + 
2(/(?/fe + ^) + 2/(j/fc -I- -I- fivk + ks)) where ki, k 2 and k^ are the estimated local slopes. 


simple approximation in (30) can be considered: since we assumed that / is locally described 


as a first order polynomial (25), it could be possible to approximate f{yk) by yj/fc where 7 is a real 


constant number that may be estimated at each time-step. 
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Vk +1 =yk + fc/(9ld+i/2) = as + hf fyt + ^ < fM >h 
Vk+I = yk + hf (yk + ^(f(yk) + hai)) . 


( 34 ) 


(35) 


As a result, (35) gives the proposed NSFD scheme (24) according to the rule #3 


of the design methodology of NSFD scheme. In the same manner as reported in 


Prop.(2.1), depending on the value of T, the total gain of |ai| could be different, that 


is why we multiply Oi by a positive constant K that should be estimated at the initial 
instants. □ 

Since f (Vk+i) — f iVk) represents the slope that locally ’’drives” the discrete differen¬ 
tial equation (|^, we are looking for a good estimation of the predicted slope through 
/ at a point y, that is located between f{yk) and f{yk+i)- Considering smooth enough 
the function f{x) where x G [|/fc, 2/fc+i], the purpose of the operator 'D{f)n is to esti¬ 
mate the derivative of / in this particular interval in such manner that an ’’averaged” 
value of f{x) for x G [yk, yk+i] can be deduced. 


Remark 2.4. This scheme (23) is the dual form of the Euler-like scheme where 
the algebraic estimator is located on the left side of 0 and is directly ’’connected” to 
the discrete values of y. In this RK-like scheme, the algebraic estimator is located on 
the right side of ^ and is connected to the discrete values of y through f. 

Remark 2.5. The lowest degree of ’’time-filtering” has been considered in the proof of 


the Prop. \2.1\\2.3\ More investigations would clarify the impact of higher degrees on 
the precision of the solution and the stability of the proposed schemes. Moreover, this 
could be an essential factor regarding the possibility of simulating dynamical models 
with noisy signals (see Rem. 


2 . 6 ) 


2.2.2 Some properties and comparison with the E-A-NSFD scheme 

The main difference in the utilization of the algebraic operator between the two pro¬ 
posed schemes is the ’’management” of the sampled y solution. In the Euler-like 
scheme, the window T is a moving window that aims to performs the derivative es¬ 
timation while the scheme runs; the window is thus initialized only at the beginning 
of the scheme. In the RK-like scheme, the window T is a static window, that aims 
to perform the derivative estimation between two steps [yk, yk-i]', the window is thus 
initialized at each time-step and the number of evaluations of / is proportional to the 


length of the window T (especially, to compute 27). 

The Euler-like scheme may require a priori few samples of y in order to compute 
V{f)n while preserving a priori the global stability of the scheme. At the opposite. 
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the RK-like scheme requests many samples due to the fact that a random process is 
involved in the evaluation of / / estimation of the slopes. In such case, the average of 
the estimated slopes is performed naturally by the hltering property of the algebraic 
estimator. 

As presented in 


2.1.3 the same type of generalization can be applied by substi¬ 


tuting the time-step h by a function (j) that follows the rule ^2 of the NSFD scheme 
design. Moreover, more sophisticated schemes can be built e.g.: association of the 
E-A-NSFD scheme and the RK-A-NSFD scheme; inclusion of the derivative estimator 
inside high order RK scheme in order to rehne to precision of the estimated slopes... 


Remark 2.6. Taking into account the properties of filtering that provide the algebraic 
estimation framework, we assume that it could be possible to simulate dynamical mod¬ 
els that involve noisy signals, like e.g. models in electronics that sometimes can not 
be dissociated from the physical noise fEBf. 


3 Conclusion and further work 

We described in this paper, nonstandard hnite-difference schemes that use the alge¬ 
braic estimation framework in order to compute an estimate of the derivatives. We 
derived two algebraic-based nonstandard schemes. The approach used for the hrst 
scheme is similar to [M], which derives a nonstandard scheme from the Caputo deriva¬ 
tive dehnition in order to solve fractional dynamical systems, and the second scheme 
aims to extend the Runge-Kutta method. Further work include: 

• extensive tests of practical problems (e.g. using the set of practical problems 
presented in [37]): 

• a complete stability study in order to characterize the stability domains accord¬ 
ing to the stiffness of the ODE; 

• the application to higher order ODEs and utilization of the dedicated toolbox 
[2T] to systematize the proposed A-NSFD procedure; 

• the application to fractional differential equations. 
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